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A multilayer shallow model for dry granular flows with the /i(/) 
rheology: Application to granular collapse on erodible beds 

E.D. Fernandez-Nieto *, J. Garres-Dfaz f A. Mangeney A G. Narbona-Reina^ 


Abstract 

In this work we present a multilayer shallow model to approximate the Navier-Stokes equations 
with hydrostatic pressure and the /r(J)-rheology. The main advantages of this approximation are 

(i) the low cost associated with the numerical treatment of the free surface of the modelled flows, 

(ii) exact conservation of mass and (iii) the ability to compute 3D profiles of the velocities in the 
directions along and normal to the slope. The derivation of the model follows m and introduces 
a dimensional analysis based on the shallow flow hypothesis. The proposed first order multilayer 
model fully satisfies a dissipative energy equation. A comparison with an analytical solution 
with a non-constant normal profile of the downslope velocity demonstrates the accuracy of the 
numerical model. Finally, by comparing the numerical results with experimental data, we show 
that the proposed multilayer model with the //(I)-rheology reproduces qualitatively the effect of the 
erodible bed on granular flow dynamics and deposits, such as the increase of runout distance with 
increasing thickness of the erodible bed. We show that the use of a constant friction coefficient in 
the multilayer model leads to the opposite behaviour. This multilayer model captures the different 
normal profiles of the downslope velocity during the different phases of the flow (acceleration, 
stopping, etc.) including the presence of static and flowing zones within the granular column. 
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1 Introduction 

Granular flows have been widely studied in recent years because of their importance in industrial 
processes and geophysical flows such as avalanches, debris or rock avalanches, landslides, etc. In 
particular, numerical modelling of geophysical granular flows provides a unique tool for hazard 
assessment. 

The behaviour of real geophysical flows is very complex due to topography effects, heterogene¬ 
ity of the material involved, presence of fluid phases, fragmentation, etc. m- One of the major 
issues is to quantify erosion/deposition processes that play a key role in geophysical flow dynamics 
but are very difficult to measure in the field. Laboratory experiments of granular flows are very 
useful to test flow models on simple configurations where detailed measurements can be performed, 
even if some physical processes may differ between the large and small scale. These experiments 
may help in defining appropriate rheological laws to describe the behaviour of granular materials. 
Recent experiments by Mangeney et al. [25] and Farm et al. m on granular column collapse have 
quantified how the dynamics and deposits of dry granular flows change in the presence of an erodi¬ 
ble bed. They showed a significant increase of the runout distance (i.e. maximum distance reached 
by the deposit) and flow duration with increasing thickness of the erodible bed. This strong effect 
of bed entrainment was observed only for flows on slopes higher than a critical angle of about 16° 
for glass beads. The question remains as to whether this behaviour can be reproduced by granular 
flow models. 

Understanding the rheological behaviour of granular material is a major challenge. In par¬ 
ticular, a key issue is to describe the transition between flow (fluid-like) and no-flow (solid-like) 
behaviour. Granular flows have been described by viscoplastic laws and especially by the so-called 
fi(I) rheology, introduced by Jop et al. [20]. It specifies that the friction coefficient /i is variable and 
depends on the inertial number / that is related to the pressure and strain rate. Lagree et al m 
implemented it in a full Navier-Stokes solver (Gerris) by defining a viscosity from the /r(J) rheology. 
They validated the model with a 2D analytical solution and compared it to 2D discrete element 
simulations of granular collapses over horizontal rigid beds and with other rheologies. Staron et 
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al. [32] and [33] applied this model to granular flows in a silo. Using an Augmented Lagrangian 
method combined with finite element discretisation to solve the 2D full Navier-Stokes equations, 
Ionescu et al. [18j showed that this rheology reproduces quantitatively laboratory experiments of 
granular collapses over horizontal and inclined planes. By interpreting the //(/) rheology as a vis¬ 
coplastic flow with a Drucker-Prager yield stress criterion and a viscosity depending on the pressure 
and strain rate, they showed that using a constant or variable viscosity only slightly changes the 
results when simulating granular column collapses of small aspect ratio. In m, Chauchat and 
Medale implemented the //(/) rheology in a three-dimensional numerical model with a finite ele¬ 
ment method combined with the Newton-Raphson algorithm with a regularisation technique. The 
numerical model was validated by an analytical solution for a dry granular vertical-chute flow and a 
dry granular flow over an inclined plane and by laboratory experiments. Previously, Chauchat and 
Medale [9] simulated the bed-load transport problem in 2D and 3D with a two-phase model that 
considers a Drucker-Prager rheology for the granular phase. Lusso et al. [ 23] used a finite element 
method to simulate a 2D viscoplastic flow considering a Drucker-Prager yield stress criterion and a 
constant viscosity. They obtained similar results taking into account either a regularisation method 
or the Augmented Lagrangian algorithm. By comparing the simulated normal velocity profiles and 
the time change of the position of the static/flowing interface with laboratory experiments of [13] . 
they concluded that a pressure and rate-dependent viscosity can be important to study flows over 
an erodible bed. Similar conclusion is presented in [22] after comparing the normal velocity profiles 
and the position of the static/flowing interface during the stopping phase of granular flows over 
erodible beds calculated with a simplified thin-layer but not depth-averaged viscoplastic model 
with those measured in laboratory experiments. 

Because of the high computational cost of solving the full 3D Navier-Stokes equations, in par¬ 
ticular in a geophysical context, granular flows have often been simulated using depth-averaged 
shallow models. The shallow or thin-layer approximation (the thickness of the flow is assumed to be 
small compared to its downslope extension) associated with depth-averaging leads to conservation 
laws like the Saint-Venant equations. These approximations have been applied to granular flows by 
Savage and Hutter m by assuming a Coulomb friction law where the shear stress at the bottom is 
proportional to the normal stress, with a constant friction coefficient /i. However, this model does 
not reproduce the increase in runout distance observed with increasing thickness of the erodible 
bed. The analytical solution deduced in [T2] proves that this system leads to the opposite effect. 
The question is as to whether this opposite behaviour between the experiments and simulations 
is due to the thin-layer approximation and/or depth-aver aging process or to the rheological law 
implemented in the model (i.e. constant friction coefficient). 

Gray and Edwards [T6j introduced the /i(I) rheology in a depth-averaged model by adding a 
viscous term. However, in depth-averaged models, only the mean velocity over the whole thickness 
of the flow is calculated (i.e. the whole granular column is either flowing or at rest). Granular 
collapse experiments and simulations have shown on the contrary that the velocity of the grains 
near the free surface is higher than that of the grains located near the bottom. During the stopping 
phase and when erosion/deposition processes occur, static zones may develop near the bottom and 
propagate upwards. The resulting normal gradient of the downslope velocity is a significant term 
in the strain rate and therefore strongly influences the fi(I) coefficient. 

To take into account the change of the velocity field in the direction normal to the topography, 
we present here a multilayer shallow model that we have developed with the n(I) rheology. This 
model consists of subdividing the domain into several layers in the normal direction and applying 
the thin-layer approximation within each layer. As a result, a velocity is calculated for each layer, 
providing a normal velocity profile. Multilayer models were introduced by Audusse |Ti] and extended 
by Audusse et al. ]4j- A different multilayer model, which takes into account the exchange of mass 
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and momentum between the layers, has since been derived by Audusse et al 0, m and Sainte-Marie 

m- 

A new procedure to obtain a multilayer model has been introduced by Fernandez-Nieto et al. 
m- Several differences appear between this multilayer model and the ones deduced by Audusse et 
al. First, in [ 13] , the multilayer model is derived from the variational formulation of Navier-Stokes 
equations with hydrostatic pressure by considering a discontinuous profile of the solution at the 
interfaces of a vertical partition of the domain. This procedure proves that the solution of this 
multilayer model is a particular weak solution of the Navier-Stokes system. Moreover, the mass 
and momentum transfer terms at the interfaces of the normal partition are deduced from the jump 
conditions verified by the weak solutions of the Navier-Stokes system. In addition, the definition 
of the vertical velocity profile is easily obtained using the mass jump condition combined with the 
incompressibility condition. 

By comparing this model with granular flow experiments on erodible beds m , [ 13] ). we eval¬ 
uate (1) if the model with the p(I) rheology gives a reasonable approximation of the flow dynamics 
and deposits of real granular flows, (2) if it reproduces the increase in runout distance observed for 
increasing thickness of the erodible bed above a critical slope angle 9 C G [12°, 16°] and (3) how the 
multilayer approach improves the results compared to the classical depth-averaged Saint-Venant 
model (i.e. monolayer model). 

The paper is organised as follows. In Section[2]we introduce the //(/) rheology and the associated 
viscosity as well as a dimensional analysis of the 3D Navier-Stokes equations. In Section [3] we 
present the multilayer approach following m to derive a 3D multilayer model for dry granular 
flows up to first order when considering the thin-layer or shallow approximation. The final g(I) 
rheology Multilayer Shallow Model (MSM) is deduced in Section [4| In Section [5] we validate our 
model using the 2D analytical solution presented in [21] and compare our results with those of 
laboratory experiments done by Mangeney et al. [23] • We show that the p(I) rheology can reproduce 
qualitatively the increase in runout distance of granular flows over erodible beds as opposed to the 
constant friction model and that the multilayer approach significantly improves results compared 
to the monolayer (i.e. Saint-Venant) model. 


2 The 3D initial system 

We consider the space variables (x, z) G fix R + C M 3 , where x = (xi, X2) G D C M 2 corresponds to 
the horizontal and z G M + to the vertical variable, the velocity u G M 3 with horizontal and vertical 
components (uh,w), the density p G M that is assumed to be known and gravity g G M. We set 
V = (<9 X1 , d X2 . d z ), the usual differential operator in the space variables, and V x := (d xl ,d X2 ), the 
reduced operator to the horizontal variable. 

The 3-dimensional Navier-Stokes equations are written as 


where g = (0, —g) G M 3 and 


V • u = 0, 

pdt.u + pS7 ■ {u® u) — V • X = pg, 


( 1 ) 


£ = —pi + T 


( 2 ) 
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is the stress tensor, with p£M the pressure, I is the identity tensor and T the deviatoric tensor 
given by 

T = r ] D(u)- T=(^ H ?K), 

\-l-XZ ± zz J 

where ijeK denotes the viscosity and D{u ) is the strain rate tensor 


where Dh{uh) 


D{u ) = Vu + (Vu)' = 


D h (u h ) d z u H + {V x w) ,s 
K (d z u H y + v x u; 2d z w 

V x uh + With these definitions, system (JTj) can be developed as 


(3) 


V. x • u H + d z w = 0, 

pd t u H + ps7 x ■ {uh <8> u H ) + pd z ( u H w) + V x p = V x • ( pD H {u H )^j + d z (j](d z u H + (V x rc)'^, 
pd t w + pu H V x w + pwd z w + d z p + pg = V x ■ (ri^(d z u H )' + V x w^ + 2d z [pd z w^j. 

(4) 

In the following subsection, the rheology and boundary conditions are presented. In subsection 
2.2, a dimensional analysis of the system is performed. 


2.1 Closures 

2.1.1 Rheology 

We consider the so-called p(I) rheology (see [20]), which is defined by 

= jA£yp_ 

1 pwir 


(5) 


where ||Z)|| = \J(D : D)/ 2, the usual second invariant of a tensor D. The friction coefficient p{I ) 
depends on the inertial number 

_ d s \\D{u )|| 

VpIps ’ ( j 

where d s is the particle diameter and p s the particle density. The solid volume fraction, denoted 
by ip s , is assumed to be constant, leading to an apparent flow density 


p = (f s Ps- 


(7) 


The variable friction coefficient is written 

/ n . P2~ Ps T 

P{I) = Ps+ T J I, 

m + 1 

where Jo is a constant value and p 2 > ps are constant parameters. Note that when the shear rate 
is equal to zero, p{I) is reduced to p s and, for high values of /, converges to p 2 - 

The p{I ) rheology includes a Drucker-Prager plasticity criterion, that is, the material flows 
when 

||T|| > p(I)p. 
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Note that the /i(/) rheology can equivalently be written as a decomposition of the deviatoric stress 
in a sum of a plastic term and a rate-dependent viscous term (see HB3): 


T = 


PsP 


J D 

l T ll < PsP 
(P>2 ~ Ps)P 


D + fjD if D 7 ^ 0, 
if D = 0; 

. Here we investigate the rheology defined by a 


with a viscosity defined as fj = . , _ 

£V5VI + IIDii 

variable friction /x(J) and a constant friction p s . Note that assuming p = p s , i.e. fj = 0, is different 
than taking fj = K , with K a non-zero constant as in [lHJ - 

The model that considers a viscosity defined by ([5]) presents a discontinuity when |jZ)('u)|| = 0. 
To avoid this singularity there are several ways to proceed. One of them is to apply a duality 
method, such as Augmented Lagrangian methods m or Bermudez-Moreno algorithm [7|. Another 
option is to use a regularisation of D(u), which is cheaper computationally, however it does not 
give an exact solution, contrary to duality methods. 

In this work, we take into consideration two kinds of regularisations of D(u). First, we use the 
regularisation proposed in [2Tj, which consist in bounding the viscosity by r)M = 250 p\fgYfi Pa-s, 
considering instead of ([5]), 

„ = - , ^ (J)P , , , . ( 8 ) 


max 


u 


/A0p\ ’ 

v M ) 


In this way, we obtain 77 = tjm if ||A ) ('^)|| is close to zero. We used this regularisation in the 
simulation of the granular flow experiments. However, as explained in Section 5.1, we cannot 
consider this regularisation in the simulation of the analytical solution, for which we take into 
account the regularisation introduced in [6j, 


n = 


M -0 p 


VI|7>WII 2 + <s 2 ’ 


where S > 0 is a small parameter. 


2.1.2 Boundary and kinematic conditions 

• At the free surface Zb + h, we consider the usual kinematic condition 

dth + u ■ n h = 0 , 


(9) 


with n h = {S/ X {zb + h), — l)/^/l + |V x (zf, + h)\ 2 the downward unit normal vector to the free 
surface. We also assume a normal stress balance 


P = PS, 

with pg the surface pressure. 

• At the bottom z = Zb we consider the no penetration condition 


( 10 ) 


un b = 0 , 


( 11 ) 
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where n b = (V x Zb, — l)/y 1 + |V x Zb| 2 is the downward unit normal vector to the bottom. 
We also consider a Coulomb type fiction law involving the variable friction coefficient p(I): 


S n b — £ n 


( K J )P] 


n n = 


U H 

U H \ 


( 12 ) 


2.2 Dimensional analysis 

In this subsection we carry out a dimensional analysis of the system i> -([12]). We consider a 
shallow domain by assuming that the ratio e = H/L between the characteristic height H and 
the characteristic length L is small. We define the dimensionless variables, denoted with the tilde 
symbol (7), as follows: 

( x, z, t ) = (Lx, Hz, (L/U)i), ( ur , w) = ( Uuh , sUw), 

h = Hh, p = p 0 p, p = p 0 U 2 p, p s = PoU 2 p~ S , 

V = PoLUfj, rj M = PoLUrj~M- 


Let us also denote 

D £ (u) = 

and the Froude number 


D h (u h ) \d z u H + e(V x u)) A 
(■ d z u H y + eVjW 2d z w 

U 


Fr = 


VgH' 


(13) 


Then, the system of equations Q can be re-written using this change of variables as (tildes 
have been dropped for simplicity): 


V x -u H + d z w = 0, 

pd t u H + pS7 x ■ (uh <8> u H ) + pd z (u H w) + V x p = V x • ( pD H (u H )) + d z [r](^d z u H + (V x w)' )), 
e 2 (pd t w + pu H V x w + pwd z v?J + d z p + pj^ = V x ■ (d z u H )' + e 2 V x «;^ +2d z (rjd z w^j. 


We also write the boundary and kinematic conditions using dimensionless variables 

• At the free surface 

&t (z b T ^) T uh\ z=zi,+h ■ V x (z b T h) ui\ z=Zb j r h 0, p Ps- 

• At the bottom 

uh ■ y x z b = w ; 

~d z u H = p(I)p^~ + O (e 2 ) . 

e \u H \ 


(14) 


(15) 
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In addition, we assume an asymptotic regime in the rheology for the friction coefficient //(/), 
namely: 

//(/) = £/J?(I). 


Consequently, 


7 ] = £ 




max 


D s (u c 


M-Op 

n m 


= £T] 


(16) 


3 A multilayer approach 


z 



Zb(x) 


X 


Figure 1: Sketch of the multilayer division of the fluid domain. 


We apply the multilayer approach proposed in [H] , Using the same notation, we denote the 
fluid domain klp(t) and its projection /p(t) on the horizontal plane, for a positive t E [0, T], i.e. 


IfA) = jx E M 2 ; (x, z ) E j. 


This approach considers a vertical partition of the domain in N E N* layers with preset thick¬ 
nesses h Q (t,x) (see figure [I). Note that X«=i = h. These layers are separated by N + 1 inter¬ 
faces T a+1 / 2 (f), which are described by the equations z = z a+ i/ 2 (C x) for a = 0,1, .., IV, x E Ip(t), 
where Zb = Zi/ 2 and z s = z^+i /2 are the bottom T^ and free surface T s respectively. We assume 
that these interfaces are smooth enough. Note that z a+ 1/ 2 = Zb + XJ 3 L 1 fys, for a = 1,..., N and 


h a = Z a+ i - Z a _i. 


Denoting £l a (t) the subdomain between T a _i/ 2 and T Q+1 / 2 and Q a (t) the lateral vertical bound- 
















ary, for a positive t G [0, T\, we obtain 

= j(x,z); x G ^f(^) and 2 a _i < 2 ; < 2 a+ ij, 
dtt a (t) = r Q _i(t) U r a+ i(t) U @ a (t), with 

@a(t) = | (x,z); x G 3/f(£) and z a _ 1 < z < z a+ i j. 


Remark 1. We need to introduce a specific notation: 

1. For two tensors a and b of sizes (n,m) and ( n,p ) , we denote by (a; b) the concatenation of 
a and b, which is a tensor of size (n,m + p). 

2. For a function f and for a = 0,1, N, we set 

4+1 ^ln a (t))|r a+ i(t) an d 4+1 ^l«a+x(t)^lr a+ i(4)' 

Note that if the function f is continuous, 

4 +§ •“ 4 ,((] ~ 4+1 — 4 + 1 ' 


3. For a given time t, we denote 

(d t z r , A 1 . V r z. 


n T,a+ 1/2 — 


a+|> V x Z a+F 


-1 


and n 


1 + 


V x z„, 1 


+ 


a+1/2 


V~Z. 


z z a+l’ 


-1 


1 + 


^c^a+1 


the space-time unit normal vector and the space unit normal vector to the interface r a+1 / 2 (i) 
outward to the layer Q a (t) for a = 0, N. 

For convenience, we write the set of equations (14) in matricial notation before applying the 
multilayer approach. First, we focus on the equations of momentum. We multiply the horizontal 
momentum equation by e, which gives 

epd t u H + ep^x ■ (u H ® u H ) + £pd z ( u H w ) + eV x p = eV x • ( rjD H {u H )) + d z (r](^d z u H + £ (V x w)' 

£ 2 {j)d t w + pu H V x w + pwdzw'j + d z p + Pj^j = V x • (S-wh) 7 + e 2 V x w;)) +2 d z (rjd z u^. 

Note that the terms involving the stress tensor (without divergence operator) are: 


eD H (u H ) \d z u H + £(y x w)^ 


V 


( d z u H y + £ 2 N x w 


2 d z w 


= V D e (u) 


£l'l O' 


0 1. 


where D e (u) is defined by (13) and I 2 is the two-dimensional identity matrix. We introduce the 
notation: 

u £ = (u h ,£w)', 


f= 0 , 


£ = 



. , 1 gH 

Wlth F+ = W 
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With this notation, we can write the momentum equation as follows 


epd t u e + epV • (u e <g> u) + V • (p£) + pf = V • (rjD e (u)£), 


and we obtain the set of equations (14) in matricial notation: 


V • u = 0, 


pd t u e + pV ■ (u e <g> u) - • (££) = —pf, 


( 17 ) 


where now X = —pi + r/D £ (u). 


In subsections 3.1 and |3.2| we define the weak solutions for our system and the process to cal¬ 


culate the vertical velocities is presented in subsection 3.3 


3.1 Weak solution with discontinuities 

Following m, we look for a weak solution ( u,p,p ) of @). We assume that the velocity u , the 
pressure p and the density p are smooth in each Q. a (t) but may be discontinuous across the interfaces 
r a+ i/2 for a = 1, ...,N — 1. Then the following conditions must hold: 

(i) ( u,p,p ) is a standard weak solution of © in each layer £l a (t). 


(ii) (u,p,p) satisfies the normal flux jump conditions at r Q+ i(t), for a = 0,..., N for the mass 
and momentum laws: 

• Mass conservation law, 

l(p\ pti) ]|r ib ^T,a+l/ 2 = 0 , (18) 


Momentum conservation law, 


(pu e - pu e ®u- -££) • n TQ+1/2 = 0, 

e J lr i(t) 

“+2 

where [(a; 6)]. denotes the jump of the pair (a; b) across T a+ i (t), 


(19) 


[ (o; 6) ]|r ,(*) = 6 )ln a+l(t ) - («; b )\n a(t) ) 


lr i (t) 
“+2 


We consider a particular family of velocity functions by assuming that the thickness of each 
layer is small enough to make the horizontal velocities independent of the vertical variable z. From 
this and the incompressibility condition in each layer, we obtain that vertical velocities are linear 
in z and may be discontinuous, that is 

W ln a (t) := U(X := i u H,aiW a ) 

where UH,a and w a are the horizontal and vertical velocities, respectively, on layer a, and the 
particular family satisfies 
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( 20 ) 


d z u H ,a = 0; d z w a = d a (t, x) 
for some smooth function d a (t,x). Note that 

= ™H, a+ i(^ x ) = u Hj a(t,x). 


Let us denote 


G Ll = d t Z a+ i + “H,q+ 1 • V*Z a+ i - W+ +i , 


G q+ i - ' V *~a+f ^ Q+ i' 


( 21 ) 


Then u satisfies the jump conditions for the mass conservation law (18) if both coincide. In this 
case we set 


<7,1 :=(?“!= G + 1} 

“+2 “+3 

Note that G q+ i is the normal mass flux at the interface T Q+ i(f). 


( 22 ) 


Moreover, using (22), the momentum conservation jump conditions (19) can be written in terms 

(pu £ \ pu E <g> u) 


of the normal mass flux as 

, ^ xZ a+ E ’ ~ 1) = 


E J|r »+j<'> 


J l r ^i C‘) 
2 


a+£» V ** a +i, X ) 


Also, by (22), 


(pft £ ; pit e <8) it) 

_ 

/ pifff 

pu H <8> uh 

pu H w\ 


lr .+i (t) 

yepw 

epuHW 

epw 2 y 


Ir . i(t) 


(23) 


12 0 
0 £ 


f pu H pu H 0 Uh pu H w\ 
Y pw punw pw 2 y 


lr «+l w 


(pit; pu 0 u) 


J l r „ , l (*) 

+ 3 


with B = e£ 1 . Using (|22|) we deduce 


(pit; pu0u) (d t z a+ i,V x z a+ i,-l ) = pG 

L J Ir i (tr 2 2 

“+3 


“+2 W|r i(t ) 


Therefore, 

(pit e ; pu E 0 u) 


= pG Q+ i B[it] 


lr ~u w 


Ir , i(t) ’ 

“+3 


a+ i, Vr^^l, 1 ) — 


\ = B (pit; pu0u) (dtz, 

> L J lr„ J .iW V 




Finally, we obtain, from the previous equality and from (23), the momentum jump condition 


^ £ \ , = 
lr a+ iw 


PG, 


a+J 


= ju it| 

2 L J lr x i«) 


1 + 


N-r ^,, 


(24) 
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3.2 Stress tensor approximation 

For a = 1 ,N — 1, the total stress is written 


, = -p . i I + T± , . 

a+| F “+2 a+|' 

ct 


where p a+ i = P + ,i = P <i is th e pressure and T* . are approximations of rj D e (u a ) at r a+ i. 


By the momentum jump condition (19), rewritten as (24), T ± , must satisfy 


1 (^ - 5T., ) £ n Q+i = t (T+ 1 -T:^ 1 ) n„ + i = 


pG <*+\ 




1 + 


V x £ 




Moreover, by consistency, we consider the following condition 


- (T + ! + T ! 
2 \ “+2 “+2 


= T 


«+§’ 


where 


T , i = r nD i i 
q+ 2 ' e,a+2 

is an approximation of ?]Zl e ( , u Q )|p . Concretely, we set 

I rv 4- i 


/ 




D h 


utr , 1 + Urj ,1 


^e,a+^,xz 


where, 


^e,a+^,xz - I 


w , i +rc 1 
“+2 “+2 


\ 


D 


e,o;+2 

J 




[it], . (25) 

1 J lr ,(t) v ’ 

OL-\- Ti 


(26) 


261 ') 


and (Qjp a+ i , Q v a+ i) is defined as follows. 

We approximate the second order derivatives in z using a mixed formulation because of the possible 
vertical discontinuous profile. We set an additional auxiliary unknown Q that satisfies 

Q - d z u = 0, with Q = (Q h ,Q„). 

And to approximate Q, we approximate u by u, a Pi(z) interpolation such that U\ z _i^ z ] +z i j = 
««■ Then Q a+ i = (Qj/ jQ+ i , Q v , a+ i ) is an approximation of Q(u) at T q+ i. 


Finally, we multiply equation (25) by e/2 and multiply scalarly (26) by vector £ n , i. Then, we 


get a linear system whose unknowns are <5 ra Q+ i. As a result we obtain 


T ± , £ ft , i = T ni ,i£ ft ,i ± - 

a+5 “+2 “+2 "+2 2 




a+i 




1 + 




lr »+“ l ’ 


(27) 


where T a+ i is defined by (|26f[)-(|26|'|) 
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3.3 Vertical velocity 


Let us recall the velocity structure requirements set in equations (20). This makes the vertical 
velocity linear in 2 in each layer. Concretely, if u a is a solution of system ( fl7| ) in n a (t), the 
vertical velocity can be recovered by integrating the continuity equation between z a _i and z E 

(z a _i,z a+ i), 

w a (t,x,z) = w+i(t,x) - (z - z a _i)v x ■ u H ,a(t, x ), for a = 1, 


Moreover, from conditions (22) at the interfaces, we obtain the relation 


«£ + i = - U H ,a) ■ VzZ a+ i + 1- 


(28) 


We therefore use the horizontal velocities deduced from the model to compute the vertical veloci¬ 
ties in the layers following the algorithm: 


From the mass transfer G 1 / 2 , which is given as data, we obtain w~[ using condition (22) at 

2 

the bottom, 

ret = uh, 1 • V x zb + d t ZB ~ Gi. 

2 2 

Then, for a = 1, ...,N and z E (z a _i,z a+ i), we set 

W a (t,X,z) = W+ 1 (t,x) - (z - Z a _l)Vx ■ U H ,a{t,x), 


w 


a+i 


= [u H ,a +1 - U H ,a) ■ V x Z a+ l + W 


where 


w 1 = w a \ = w + - h a X7 x ■ u H a . 

OL+ 2 ir a+l/ 2 ^) a—o 


(29) 


In this way, the velocity vector u is the piecewise smooth function such that u(t,x,z) |n a (t) = 
u a (t, x, z) for a = 1,..., N, where 


U a (t,X,z) = (uH,a(t,x), W+_i(t,x) ~ (z - 2 a _ 1 )V X • U H ,a(t, x )^j , 


(30) 


and w + _ 1 (t,x) is computed using (29). 


4 Weak solution of the first order model 

In this section we derive the model of order e. 
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4.1 Pressure 


Using the vertical momentum equation in (|14|) we can deduce an expression for the pressure. We 
write this equation up to order e 2 for each 


ayer: 


1 


<9 zPa = ~Pjr 2 + • ( rjdzUH,a ) + d z (2 T]d z w a ). 

Taking into account the requirements d z UH,a = 0 (see (20)) and 77 = erf (see ©), we get: 

d z p a = -/°^2 + £dz ( 2l fdzWa) • 

Then, we obtain the hydrostatic pressure framework (up to order e) in each layer: 

1 


d z p Q = -p 


Fr 2 ' 


Now, by the continuity of the dynamic pressure (see S3!), we can deduce that 

Pa(z) =Ps + (z b + h- z), 
where ps is the pressure at the free surface. 


(31) 


4.2 A particular weak solution 


Noting that u a is a weak solution of the system (17) in Q a (t), let us consider the weak formulation 
of (17) in fi a (t) for a = 1, Assuming u a £ Lr( 0, T ; i7 1 (0 Q ,(t)) 3 ), dt.u a € L 2 ( 0, T ; L 2 (Q a (t)f) 
and p a £ L 2 (0, T; L 2 (Q a (t))), then a weak solution in Q a (t) should satisfy 


4 / 

6 Jn a (t) 


0 = 


pf ■ vdfl = 


/ (V • u a ) <p d£l, 

/ pdtu £:a ■ v + / p( u e>a ■ S7u a ) ■ v dfl + 
'Q a (t) JQctit) 


(32) 


1 

£ Jn a (t) 


(v • (p a s)) -vdn--[ (v • (r a £)) ■ v dn, 

la(t) £ 

for all ip £ L 2 (Q a (t)) and for all v £ i7 1 (0 Q ,(t)) 3 . 


We consider unknowns, velocities and pressures, that satisfy (20) and the system (32) for test 
functions such that 

d z ip = 0 

and 

v(t, x , z) = (v H (t, x), (z- z b ) V ( t , x)) , U| a/j?(t) = 0, (33) 

where Vh and V ( t , x) are smooth functions that do not depend on z. 


We will now develop (32) in order to obtain the mass and momentum conservation equations 
that satisfy the weak solution for this family of test functions for each layer. 
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Mass conservation 


Let ip = <p(t,x) a scalar test function and u a a weak solution of (32), from the mass conservation 
equation we get 


0 = / (V • u a ) (p dfl = 


<p(t, x) ( j +2 (V x • uh.q + dzWa'jdz | dx = 


'i F (t) 


/ <p(t,x) I 

JI F {t) 


rz . i 
a +? 


Vx ■ UH,cedz + W J — W + ! dx, 


for all a = 1,..., N. Using Leibnitz’s rule 


r*+h 

J Z 1 
a ~7 


’ 'U j H,ot.dz — V a 


r *«+i 


' 2 1 
a-2 


I ' V ^a+i - V xZ a _i, 


~ H , a +± ’ x ~«+2 ' 


and this leads to 


/ 

J i F (t) 


<p{t, x) ( V x • ( h a U H)0l ) - Utf )Q: • V X Z Q+ 1 + W Q+i + ' V X Z Q _1 - ) dx = 0. 


Taking into account (|22|) and <9fh Q = 9*2 : q+ i — dtz a _i we obtain the equation 


/ <p(t,x) d t h Q + V x • (h Q u H , a ) ~ G a+ i + G q _i dx = 0, 

d/ F (t) \ 2 2j 

for all y?(t,.) G L 2 (Ip(t)). Thus we get the mass conservation law for each layer 


dth a + Vi • (h a UH,a) — G ,i — G a _ i , a — 1,..., N 


(34) 


where £^+ 1/2 and G 1/2 stand for the mass exchange with the free surface and the bottom respec¬ 
tively and both should be given data. 


Momentum conservation 

First we develop the variational formulation of momentum equation taking into account that 
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V • ( T a £) ■ v (in = - [ (T a ■ S) : Vv dSl - 

J VL a (t) 


,r - + lW 


T , X E ) v ) • ft a+ idT + 


Now we can write the momentum equation as 
1 / 




pf-vdQ = 


£ Jn a (t) ^ 

+ 1 / 

£jr „aW 


iW 


/r lW vv T ^)")-^ dr - 

Q_ 2 V 




: V# dQ + - f 

^ JQ a (t) 

-P a+ i_£ + T a+l £ 


i+ 

a-1 


• Vu a J 

• v dQ — 

(T a £) : 

: Vu dQ + 

"«+i) 

■ v dr- 

™--i) 

■ v dr. 


(35) 


Let u e H l (rt a ) be a test function satisfying (33). We develop the momentum equation in 


(35) by integrating with respect to the variable z and by identifying the horizontal and vertical 


component of the vector test function v. In addition, taking into account the hydrostatic pres¬ 
sure framework, we can leave out the equation corresponding to the vertical component. This 
is equivalent to considering the vector test function where the vertical component vanishes, i.e. 
v = 0). Therefore, the horizontal momentum equation reads, for a weak solution u and 


for all a = 1,..., IV: 


— [ pf ■ {v H ,0)dQ= [ pdt(uH,a,£w a ) ■ {v H ,0)dn + 

£ JJ n a (t) 

+ / p[{u H)a ,£w a ) ■ V(u H ,w a )) ■ (v H ,0)dtt - 
Jn a {t) v 7 

- 1 / (p a £) : X7(v H ,0)dn+ 1 [ (T a £) :V(v H ,0)dn + 

£ Jn a (t) £ 


+ -/ 

6 ir „ o(t) 


(36) 


~P a +k £ + T a+ l £ ) n a+\ ) ' (*H, 0) dT - 


e J r i (t) 


~P a -h £ + T ^± £ ' {v H ,0)dF. 


We develop each term of this equation, taking into account that 

d z u H ,a = d z v H = v\ = 0. 
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z a + l/1 


/ pd t UH,a • vh d,Cl = / / pd t u H ,a ■ v H dzdx = / ph a d t U H ,a ■ Vh dx. 

Jn a (t) di F (t) J J ip (t) 

z c,-l/2 

/ PyUHyCt ‘ ^xP^H.aj ' V H d£l — / p (h a UH,o ' V x.UH.a) ' V}jdx. 

Jn a (t) v 7 


/ z a+l/2 


iQ. a (t) 


Pa Vj ■V H dQ = 


Ip(t ) 

/ z c* + l/2 


dz V x • vh dx = 


\fa-l/2 


V., 




/ z a + l/2 


p a dz ■ vndx + 


\ z a-l/2 
/ z a+l/2 


'9/fW 


Pa dz vh • ndr = 






_ , dz 

V xPadz + Pa~T 
dx 


Z . 1 
“+2 


\^q:-1/2 

vndx = (*). 


1/2 


Moreover, 


'a + 1/2 z a + l/2 

J V xPadz = J V x (ps + -^2 + h - z)j dz = 


z a-l/2 

z a + l/2 

-! 

z a- 1/2 


z ct-l/2 


(^xPS T p 2^ x dz — haVxPs + p 2^°^ x (^6 T d) • 


Then, we continue the computation of (3). We obtain 
(*) = - J ^ ( h a VxPs + (z b + h) + p a+ iV x z a+ i - p a _iV x z a _i^j ■ v H dx = 

= - J (haVxPs + jr^haVx (z b + h)J • v H dx - 


/ P a+ in a+ y(v H ,0)Jl + 
Jlppt) 2 2 V 


V 1 + 

^I Z Q+i 

dx+ [ p a _in a _i ■ (v H ,0)\ 1 + 

V 


d/ F (i) 2 2 V 


V x z n 


'Ha(*) 


/ z a+l/2 \ 

T H.a '■ V X V h dQ= / Th, a ■ ^ xV h dz \ dx = 

dip (t) J ) 

\ z a -l/2 / 


2 

1 dx. 

2 
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11 F {t) 


/ z a + l/2 \ / Zq + 1/2 \ 

/ Th/ )Q dz : V x v H dx = - / V x • / Ttf jQ dz • v H dx + 

J I d/ F (i) J I 

\fa-l/2 / \fc«-l/2 / 


\ 2 a- 

/ 2 a + l/2 \ 

+ [ [ T H , a dz J • • ndT = - [ 

J dip if) J I dIp(t) 

\ol- 1/2 / 


/ z a + l/2 


v x • / dz • v H dx. 

” / JI F{t) J 

\a-l/2 / Vfa —1/2 

Because d z un, a = 0, we obtain that ||-D e (it Q )|| is independent of z up to order e, since 




/ D H (u H ,a ) 0 

\ o 29, W , 


Then, from (16), we get 

/ ' T H) a cfe = 

J Z 1 


+3 £r] 0 D H (u H}Ol ) dz = O(e). 


Therefore, we can neglect the term V a 


r«+h 

J Z 1 


T H, a dz , which corresponds to the horizontal 


diffusion, since we are interested in the first order model. 


/r Q+ iW 


' 1 Fit) 


-p a+ i£ + T a+i £ ) n a+ 1 ) • (t? H ,0) dT = 


-p a+ i£ + T a+i £ ) n a+ 1 ) ■ 0)^/1 + 




dx. 


Introducing these calculations in (36) and taking into account that / = (0, , we obtain 

(ph a dtU H ,a T ph a UH,a ' ^ x^H, a T h a \7 xPS T «Vj; (z fe + /l)) -V H dx + 


'IFit) 


+ 4 w l p ” + ^' ( ”"' 0) ' /I + 


^x Z a+\ 


~P a - in a _i • (vh,0)\ 1 + 




dx + 


+ 1 / 

£ JI Fit) 


-p a+ i£ + T a+ ,£^j ft a+ -(v H ,0)^l + 


V x z. 


- ( ( if + T+if ) n Q _i ) • (ifa,0Vl + 


Q 2 Q— 2 ) 01 2 


V -T z,.. 




dx = 0. 
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Note that -p a+ i£n a+ i ■ {v H ,0) = p a+ in a+ i ■ (v H ,0), then 


J (ph a dtUH,a + phaUH,a-'V x UH,a + haVxPS + J^-2h a Vx(Zb + h)')-VHdx + 


+ 


’ lF{t) 


- T a+ ^™«+\ ) ' (VH,0)\/1 + 


^* Z a+\ 




VxZ, 


x z n -l 


dx = 0. 


Moreover, 


T a+ i £fi a+h ) ' (vh, 0)^/1 + 




eT jt ! T ! 

H,a+ 2 xz,a+2 


eT 


T~ 


^x z a+ l 

-1 


V H 

0 


= (^,a + I V ^ + i- T ^ Q+ i ) 


Therefore, 


//_p(t) L 


ph a d t UH,a + ph a UH,a • ^xU H , a + h a X7 x p s + (zb + h) + 

vndx = 0. 


+ ( T if, a +i V *^+| e T ^,a+i) £ T ia-| 


And this yields, for each layer a = 1, ■ ■ ■ , At, the momentum equation 

ph a dtUH,a T ph a U}{ a • V x^H,a T haS7 xPS T J7^2 ^o(N x ( z b T h) T 


+ T 


H,a+\^ x ~ a +\ £ i iz,a+t 


- -T“ 


1. 


- ( T+ , V x z„ i - -T + , 

H,a— | a— 2 ^ xz,a— | 


= o. 


Remark 2. Observe that 


T , V r 2 , i-T , 

H,a+| x «+5 e a; 2 ,a+| 


-T Tn 


e “+ 


i crt o+i 


1 + 


N ;r Z 


a+i 


J 77 


where [ ■ ]h denotes the first component. Now by (21) we obtain 


1 

/ 


2 


1 ~ / 


-'J’ il’n i i \ 
e "+2 " + 2 1 

h 



H 

- r „++"«+iV 1 + 

V *~a+i 


2 1 

-pG„i full 

2 a+ 2 1 lr a+ i(t) 


= T H,a+ ± V ^^o+i - \ T xz ,a+i “ ^ G a+i (**#.«+! “ S H>«) ’ 


2 2 
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and analogously 


-T+ i £n i 

rv-- o 


1 + 




= T H,a-^xZ a _i - ^T xza _ i + ^pG Q _i (Ufl,a - 


J H 


Remark [2] allows us to re-write the momentum equation 

ph a dtUH,a “k ph'a'U'H,a ' ^x^H,a ~k ha^xPS ~k aV x (^6 + h) + 


+ I r H,a+i V:cZ a+i £ ^,a+i j f T H,a.-\ £ T xz,a-\ 1 “ 


= ^pG a+ i (wj?,a+l - "Rff,a) + ^pG a _i (UH,a - UH,a-l) 


By combining the previous equation with (34) we get 


pdt ( h a UH,a ) “k P^x ' (h'a'U j H,a ® UH,a) ~k ha^xPS ~k Ft 1 k aV x (zj, + h) + 


+ I ^Va+4 V ^ Q+ i ^,a+| 


— — iVx2.— I-TL 


tf,a-± v x ^a— 


2 e 


= ^pGq,_)_I (m/Z>-|- 1 + W/f,a) - ^pG a _ i (Wff, a + UH,a-l ) 


Note that 

f 


H,a+^^ x ^a+^ £ '^'xz,a -|-| 


-T Q+ ifn Q+ i A /l + 




J H 


V°D a+ i£n a+ i 


-1 / 


.W 1+ 

^ Z a+§ 


o n I “#■“+* + W R«+| , ^ 

= <+ h DH ( - 2 - I '^+1 


- v. 


it; +iy 


1 T CO 1 
4 ™_j_ 1 


-ri a+ lQ H ,a+ 


1 — 
2 


+0(e). 


We define 


K 


a+| = ~zV nJ .iQ 


where ry° + i is a first order approximation of 77 at z = z a+ 1 . We obtain 


a+i^H,a+i> 
a+1 


(37) 
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max ( ||Q 




MJ„+i)p a+ i 


Vm 


(38) 


with 


N 

P a+ l 2 =PS+Yr 2 V 

/5=a+l 


/ , i = 


^IIQ^a+i I 

V /p «+i /ft 


(39) 


These expressions of r? Q+ i and / a+ i are obtained from definitions (J8J) and (|6j), respectively, by 

considering the hydrostatic pressure approximation (31) with the definition of p Q and with the 
following first order approximation of ||Z)(w)|| at z = z a+ i, 


\\D(U)\\ U=Z , ~ \\Qh,o.+ \\V 


(40) 


By re-writing the momentum equation again, we obtain up to order e, 

pdt x ' ® 'U j H,ol) H" ha^xPS “1“ x (Zb T h) — 

= K a -\ - K a +\ + 2 p( ^ a +h + UH,a) - ~^pG a _i ( U H ,a + UH,a- 1 ) , 

the horizontal momentum conservation laws, for a = 1, N. 


(41) 


Remark 3. We must impose friction at the bottom (T i with a = 1). We can translate (15) 
into the notation of the multilayer approach, giving 


w In = 0 ’ 
2 


~Vi Q h ± = 


u 


H 4 


£ 2 


U H,\ 


Therefore to impose the friction condition, we should change definition (31) of Ki, taking into 

I I 2 

account that 

h +1 = u H , 1 - 

Then we obtain 

_» 1 „ _» Ti.TT 1 

(42) 


I<i = - 1 r]°iQ H i = -fi(I)°pi 1111,1 


£ 2 


\ u H,l\ 
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4.3 Final Model 


We have obtained the dimensionless final system given by (21), (22), (|34|) and (|3T|)- (42). The last 


step is to return to the original variables taking into account subsection (2.2). We obtain the final 
multilayer system, for a = 1,..., N, 


dtK T Vx • ( h a iiH,a ) — G q+ i G a _i , 

pdt (h a UH,a) T Vj, • {ph a llH ta ® UR t0 ^ T 

+ h a \7 x p s + pgh Q V x (z b + h) = K a _ i - K a+ 1 + 

+ -pG Q+ i (UH,a +1 + Wif.a) ~ ~^pG a _i (UH,a + Ur, a- 1) , 


(43) 


where 


G «+i = 3t*a+± + ■ VxZ a+ l - W+ + i = 3t* a+ i + U H , a ■ V x Z a+ 1 - W q+ i • 


and 


X a+i = "VlQ 


2 '^5 Q! +2 




2JH 


1 + 


^x Z a+\ 


(44) 


Ki = -g{I)pgh ™ H;i 
2 \u H . i 


System (43) must be closed by setting the vertical partition of the domain. For this, we can 

h where 


write the thickness of the pre-set layer based on the total height. That is, we set h a = l , 
l a is a positive constant, for a = 1, • • • ,N, and 


N 


Ct=l 


= 1. 


Note that G a+ i can be written, by summing the mass equations from 1 to a, as 


G 


a+ 


1 

2 


Gi + ^ {d t hp + Va; • (hpURfi)). 
p=i 


(45) 


Moreover, for the special case a = N and assuming no mass transfer with the atmosphere, i.e. 
G n+ i = 0, the above equation leads to 


8th + Va; 


N 


h ^ IpUR^ 


.4=1 


Gi. 


By introducing this in (45) we obtain 


a / N 

G a+ l — G l T ^ ( Ip ( Va; (hllR p ) ^ ( V x ‘ (l'yh'UH,'y') 

2 4=1 V 7=1 



(46) 
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Let us define L a := l\ + ■ ■ ■ + l a and £ a>7 = y^((fg 7 — 1/3)1 7 , where dp^ is the standard Kronecker 

/3=i 

symbol. That is, 


£a,7 — * 


(l - (/i H- bl a ))l 7 , if 7 < a, 

—(h + • • • + Z Q )Z 7 , otherwise, 


for a, 7 £ {1,... ,N}. Then, we can write the mass transfer (46) in the interface as 


N 


G aJ< _i — (f L a ) G l + ^ ^0,7^a; • (/rti// j7 ), 


(47) 


7=1 


for a = 1,... ,N. 


Next, considering q Ha = hun,a and after some straightforward calculations, the final system 


(43)-(44) is re-written, for a = 1, • • • , N, as 
dth + V x • l/3q H p \ = — Gi. 


L3=i 


dtqH,a + V a 


Qh.cx ® Qh,c 

h 


+v„i4+—) - ™ dxh+ 

2 p J p 


+ y~] -jJJ- f {o.H,a + QH,a- 1) Ca-1,7 ~ (<h/,cH-l + Qh,o) Cayyj ^ ‘ (gff, 7 ) — 

= ^ (9//,cH-1 + QH,a) (1 - L a ) — (q Ha + qn,a- 1 ) (1 ~ L a - 1 ) j Gl — 

^ 0 


-ghV x z b + — ( K a _ i - K a+ 1 


4.4 Energy associated with the final model 


In this section we study the energy balance of the obtained system (43)-(44). 


Theorem 1. Denoting the energy of the layer a = 1, • • • , N for the system (|43|)-(|44|) by 

E a — h a 


u H,a\ .PS . . 

+ — + 9 1 z b + 2 


X 2 p 

the following dissipative energy inequality is satisfied: 

h 


N 


pdt I E a J + pV : , 




N 


.a=l 
2 


y. u H ,a ( E a + pgho 


< h a d t ( PS + pgz b ) - 


r I -. I / t\ \UH,N\ 

-pgh \u H) i p(I) --- Vn+- ~ 

h N 


TV—1 

E 

a=l 


{l^H,a+l H a ) 


v a+ i - Gi (ps + pg(zb + h)). 
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Proof.- 

Firstly we multiply the momentum equation for each layer by UH,a and use the mass equation 
to simplify the convective terms. Secondly we sum up the obtained equation for layers a = 1 to 
a = N and then we obtain that the global system has a dissipative energy balance. 

Now, we write the momentum conservation equation in terms of the velocity using the continuity 
equation 


ph a d t u Hi a + P (V x u H ,a) h a u H ^ a + h a v x p s + pghoSx Oft + h) = 

= -fr'a-i ~ + nPO+- (uH,a+\ ~ Ur, a) + 7,pG a _ I {UH,a ~ UH,a- 1 ) • 


(48) 


Multiplying (48) by UH,a we obtain 


ph Q U H ,a • d t UH,a + P {(^xU H ,a) h a U H ,a ) • U H ,a + h a U H , a ' V x (p S + pg Oft + h)) = 


= U Hy a ■ ( K a _ l - K a+ l J + -pG a+ l ( UH,a+ 1 ' Uff.a ~ \uH,a\ ) + 


+ ^pG a _ 1 (\u H)a \ 2 - U H ,a ‘ UH.a-i ) • 


This gives 


h'a'U'H,a) ' 'U j H,a — ^a 


\UH,a | 


ha^H,a — 


= v x 

which can be re-written as 


\tiH,a\ 2 J - \ lu^al 2 , - , 

il a UH,a I rj V x \h / a / U'H,a) ? 


/i^ i2 \ i2 

ph a UH,a * dt'U>H,oi. P^x * ( ^ h'Oi'U'H,a | P ^ V# * {h a UH,a) H - 


+h a U H ,a ■ V x (ps + pg Oft + 0) = ' ( -Kq,_ I ~ K a +±- ) + 


(49) 


+ (‘UH,a +1 • I (|itff,a| 2 ~ ' Uh,ol- 1 


Mg, 


Let us consider the mass conservation equation multiplied by p — -1- ps + pg(z b + h) , 


1 ^ |2 |12 
'll TT 'll TT 

p- — Y^-d t h a + p —Vj. • (. h a u H , a ) + (ps + p^Ob + L)) d t h a + (ps + pg(z b + /i)) V x • ( h a u H)0i ) 


= P 


\UH,c 


(g q+ i - G a _ij + Os + pg(z b + h)) (g q+ 1 - G a -§) . 


(50) 


24 












Summing (49) and (50), we obtain the equation 


pd t 


| uh, ( 


" ] + p\7 j 


| U H ,c 


-h a u H , Q + (ps + pg(zb + h))d t h a + 


+ Va; • ( h a (ps + pg (z b + h )) U H ,a) = U H)0l ■ (K a _ 1 - ^a+\) + 


+ G a+ i ( P 


+ ps + pg(z b + h)) - G a _ i (p 


+ ps + pg{zb + h) 


We can reformulate this as 


pd t 

+ pV 


hr 


!%L + r? +s(2l + fc) j 


- - pgh a d t z b - pgh a d t h + 


ha +^j + g( z b + h) \ U H ,a 


= u a ( K a _ i - K a+ 1 ) + 


+ G a+ 1 UH,a +ps + pg{z b + ^ 


G q _ i I p UH ’ a ™ H ’ a 1 + ps + pp(z b + h) 


Note that 


-pgh a d t h = -pghadt ^ - pgh a d t ^ = 


-ps-di + pghdt— - pgh Q dt - , 


then, 

pd t 


u I \u H ,a\ , PS , ( , h N 

2~ + 7 +5 V* 2, 


+ y ( hd t h a - h a d t h ) - /lamps' - pgh a d t z b + 


+ pV 2 


i I \u H ,a\ PS ( , h\\ _ I h. _ 

I -^-1 - ~ + 5 ( £& + Y J I + gh a -U H ,< 


= V„(K a _i-K a+ 1) 


+ <7+ i ^p Mg,a+1 2 Ujf,a + ps + pg{z b + i ™ H ' a 1 + ps + P9( z b + h )^j 


Denoting 


Ea — h a 


l U H ,: 


2 —+ 7+«(^+^ 
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we have for a = 1 ,N the following energy equality 


pd t E a +y (hd t h a - h a d t h) + 


(i) 

+ pVa; 


( 2 ) 

E a + ) U H ,c 


= u a (K i - K i i ) +h a d t (ps + pgzb) + 


+ G 


a+\ 1 P~ 


(3) 

UH,a +1 ' UH,a 


(4) 


+ ps + pp(z& + fr) ) -G I ( p Mg,Q 1 + ps + pg{zb + M ) • 


(5) 

(51) 

Now we sum up (51) from a = 1 to a = N. We take into account that ^ Q=1 h a = h- G n+ i = 0 

(there is no transfer with the atmosphere) and uh,o = uh,n +l = 0 (velocity of the bottom and 
atmosphere respectively). This gives, term by term: 

/ N \ 


(i) P d t Y, E o 


^a=l 
/ N 


N 


(2) hdt I ha ) ~ ^2 hadth = hdth — hdth = 0. 


\Ot=l 

N 


a=l 


(3) pV a 


.a=l 


V u H ,a ( E a + pgha- 


Uh ,1 


(4) Taking into account K a+ i = —r] a+ iQ H a+ i and K i = — pgh ’ p(I) when we sum all the 

2 2 ’ 2 2 

layers, we obtain 


Af-l 


-pgh\ uh,i\ p(I) + UH,NP N+ i.Q H}N+ l - ^2 (uH,a+i-u a )v a+ iQ H ,< 


a=l 


where Q^ a+ i is an approximation of 8 z uh in T a+ i. We consider 


Q 


_ UH,a +1 — Uh )C 
H ’ a+h “ h a+ t 


Q 


H,N+h ~ 


uh, n+i - UH,N 
hN 


with h a , i being the distance between the midpoints of layers a and a + 1. This gives 


-pgh \u Ht i\p(I) - 


\Uh,N\ 

hN 


> iV-1 x2 

E \ u H,a+l ~ u a) 

h : ' q a+ i 


a= 1 


K+i 


which is a dissipative term. 


26 




















(5) Considering G N+ 1 = 0, we have 


~G i (p s + pg(z b + h)). 
2 

Finally, by summarising (l)-(5), the proof is completed. 


□ 


5 Numerical tests 

The numerical approximation is performed in 2D (downslope and normal directions). We re-write 
the model as a nonconservative hyperbolic system with source terms as in m ■ Then a splitting 
procedure is considered. First, we set aside the term that appears in the internal interfaces and a 
standard path-conservative finite volume method is applied. These path-conservative methods were 
introduced in |26j . To deal with the Coulomb friction term, we use the hydrostatic reconstruction 
introduced in [2], which is applied in [8] to solve the Saint-Venant system with Coulomb friction. 
The main advantage of this reconstruction is its great stability. 


The second step is to solve the contribution of the term in the internal interfaces, which rep¬ 
resents the mass and momentum exchange between layers. In this step, a semi-implicit scheme is 
employed, taking into account the regularisation of ||Z)(’U Q ,)|| mentioned in Section 2.1.1 in order 
to avoid the singularity when ||D(w a )|| vanishes. 


In order to validate the Multilayer Shallow Model (denoted MSM hereafter) with the fi(I) 
rheology, we compare it to (1) a 2D analytical solution for steady uniform flows over an inclined 
bed and (2) laboratory experiments of granular collapses over an inclined plane covered by an 
erodible bed made of the same material. 


5.1 Analytical solution 

Let us first compare the model to the 2D analytical solution deduced in |21| for a uniform flow 
over an inclined plane of slope 6 and thickness H > 0. This solution is obtained by imposing zero 
pressure and zero shear stress at the free surface and a no-slip condition at the bottom. 



Figure 2: Sketch of the analytical solution. 
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By denoting u and v the downslope and normal velocities, p the pressure and r the shear stress 


and by taking the rheological parameters defined in Section 2.1.1 the analytical solution reads 

Vvsgcosd (if 3 / 2 — (H — z) 3/2 ) , 


_ _2 / tan9 - p s 

3d 0 \ A p, — tanO + p s 


u(z = 0) = 0, v = 0, 

p(z) = pgcosd (H — z ), 

t(z) = p{I)p = pgsind (H — z ), 

p(z = H) = 0, t(z = H) = 0, 

„ p(I) = tan{9 ), 


for 2 E (0, H). 


(52) 


For the numerical simulation, as in the analytical solution, we consider a uniform flow with 
constant thickness H = 1 nr and velocity u = v = 0 m- s -1 at the initial time t = 0 s. The 


boundary condition at the free surface and at the bottom have been set as in (52). At the right 
and left boundary, we use open boundary conditions. 


Note that at the free surface we have 

p = 0 and ||P(„)|| = = 7 ( A „*t’ ~) = 0. 

As a result, we cannot use the regularisation ([8]) since its denominator 

max f\\D(u)\\, 

\ VM ) 

vanishes at the free surface. In this case we use the regularisation 


V 'l|PWII 2 + <5 2 ’ 

where 5 > 0 is a small parameter (see 0). 


We choose the rheological parameters /o = 0.279 and p s = 0.363 ~ tan(19.95°), p, 2 = 0.74 ~ 
tan(31.8°) and the particle diameter d s = 4 cm with solid volume fraction ip s = 0.62. The slope 
angle is taken as 9 = 0.43 rad ~ 24.64°. Figure [3] shows the good agreement between the simulated 
and exact solutions for the profiles of the velocity, pressure, shear stress, p{I) and ||Z?(it)||. It also 
shows the downslope velocity at the free surface as a function of the slope angle. These results are 
computed using 50 layers in the MSM. 


Figure [4] shows the computing time required to simulate 50 seconds and the relative error 
between the computed velocity and the exact solution using a different number of layers. The error 
is computed by 


A u 
u 


N 


2^i=o\ u ex,i 'U'sim,i . 


2^i=0 


U 


2 

ex,i 


2 


(53) 
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Figure 3: Comparison between the analytical solution (dashed and solid lines) and the simulations obtained using the 
MSM with the rheology (symbols), (a) Analytical and simulated downslope horizontal velocity u, pressure p and strain 
rate ||Z)(w)||; (b) Analytical and simulated shear stress and friction coefficient (c) Comparison between the simulated 

(symbols) and the exact (dashed line) horizontal velocity at the free surface as a function of the slope angle. 


where u ex (respectively u S i m ) is the analytical (respectively computed) velocity and M is the 
number of partitions of the mesh in the horizontal direction (in this case M = 20). Note that for 
slopes smaller than arctan(/i s ), the surface velocity is zero because the mass does not flow. The 
error decreases as the number of layers increases and is less than 10% for 20 layers. The main error 
occurs near the free surface where the gradient of the horizontal velocity is large. 

5.2 Comparison with laboratory experiments 

We will now use the multilayer shallow model to simulate the laboratory experiments performed 
in [213 • The objectives are threefold: (1) to evaluate if the model with the /i(/) rheology gives a 
reasonable approximation of the flow dynamics and deposits of real granular flows, (2) to observe if 
it reproduces the increase in runout distance observed for increasing thickness of the erodible bed 
above a critical slope angle 9 C e [12°, 16°] and (3) to show how the multilayer approach improves 
the results compared to the classical depth-averaged Saint-Venant model (i.e. monolayer model). 
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Figure 4: (a) Computing time as a function of the number of layers in the MSM and (b) relative error between the 

computed and exact velocity for simulations over a slope of 9 = 24.63°. 



Figure 5: Sketch of the initial and final state of the granular collapse. A granular column with a thickness ho = 14 cm 
and a length ro = 20 cm is released on an inclined plane of slope 8. The plane is covered by an erodible bed of thickness 
hi made of the same material. When the flow stops, the maximum final thickness is hf and its final extent rf. 


The variable rf denotes the runout distance, i.e. the length of the deposit measured from the 
position of the front of the released material at the initial time located at x = 0, tf denotes the 
flow time from t = 0 s to the time when the material stops and hf denotes the maximum final 
thickness of the deposit (see Figure [5]). 

In the laboratory experiments performed in [24] . subspherical glass beads of diameter d s = 0.7 
mm were used. They were cohesionless and highly rigid. The particle density p s = 2500 kgm~’ i 
and volume fraction p s = 0.62 were estimated, leading to an apparent flow density p = ip s ps = 
1550 kgm ~ 3 . 

In order to use the p(I) rheology, the rheological parameters (p s , p 2 and Iq ) are taken as in 
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[18] . according to the measurements made in the experiments of [23 j and [29] . where the effect of 
lateral wall friction is taken into account empirically. These parameters can be obtained by fitting 
the curve h stop (9), where h s t op is the thickness of the deposit lying on the slope when the supply is 
stopped after steady uniform flow (see [28] for more details). As a result, we take here Iq = 0.279, 
as in [20], n s = tan(25.5°) ~ 0.48 and y s = 0.74 ss tan(36.5°). 

This experiment has been simulated for different slopes 0 and thicknesses hi of the erodible 
bed: 9 = 16° and hi = 1.4, 2.5, 5 mm, 9 = 19° and hi = 1.5, 2.7, 5.3 mm, 9 = 22° and 
hi = 1.82, 3.38, 4.6 mm, 9 = 23.7° and hi = 1.5, 2.5, 5 mm. Note that the model does not take 
into account the effect of removing the gate during the initial instants even though it has a non- 
negligible impact on the flow dynamics as shown in [18]. For instance, when the gate is taken into 
account, even with no friction along it, the flow is substantially slowed down however the deposit 
is almost unchanged. All the simulations are performed using 20 layers. 

We compare hereafter (i) the constant and variable friction rheologies and (ii) the monolayer 
and multilayer approaches. 

5.2.1 Deposit profiles 

Let us compare the deposits simulated with the /u(I) rheology and with a constant friction coef¬ 
ficient n s for different slopes 9 and erodible bed thicknesses hi. Figure [6] shows that the deposit 
calculated with the variable friction coefficient y(I) is closer to the experimental deposit than the 
one calculated with a constant friction coefficient y s . The runout distance with the constant co¬ 
efficient fi s is always too long except at 9 = 19° and hi = 5.3 mm (see Figure [6jl). To properly 
reproduce the runout distance with a constant friction coefficient, we need to increase its value. 
For example, with a slope 9 = 16° and an erodible bed thickness hi = 2.5 mm (Figure [6^,), we 
need to use the value y s = tan(27.3°) to produce the runout observed in the laboratory experiments. 

Figure [7] shows, for a slope 9 = 22° and hi = 1.82 mm, the final deposit obtained using 
the constant or variable friction coefficients for both the multilayer and monolayer models. The 
difference between the multilayer and monolayer models is stronger when using the /u(J) rheology. 
For instance, the monolayer approach changes the full deposit profiles for the /x(7) rheology, while 
it only changes the front position for /i s . The multilayer approach makes it possible to obtain a 
deposit shape which is very close to the experiments with the /r(J) rheology. More generally, the 
shape of the deposit is closer to the observations with y(I) than with fi s in the Multilayer Shallow 
Model. 

5.2.2 Effect of the erodible bed 

Figure shows two zooms, one near the front (I) and one near the maximum thickness of the 
deposit (II), for 9 = 22° and different values of hi (see Figure [7] for the approximate location of 
these zooms). With the variable coefficient /u(I), the runout distance r/ increases as the thickness 
of the erodible bed hi increases (see Figure l8b(II)) as observed in laboratory experiments. On 
the other hand, with a constant friction coefficient y s (Figure [8^(11)), the runout distance iy 
decreases with increasing hi. Note that in both cases the maximum final thickness hf decreases 
with increasing hi as in the experiments (Figure |8^i(I),b(I)). 

Figure [9] shows that the decrease in runout distance with increasing hi for constant friction y s 
is observed for all slopes, e.g. 9 = 0°, 10°, 16°, 19°, 22°, 23.7°. For the /r s -model, the multi- and 
monolayer models follow the same trend. Note that this nonphysical decrease in runout distance 
with increasing hi has been demonstrated analytically in |12j for the monolayer model. Moreover, 
laboratory experiments show that when the thickness of the erodible bed increases, for slopes 
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Figure 6: Deposit obtained in the experiments (solid-circle blue line) and with the Multilayer Shallow Model using a 
constant friction coefficient p s (dotted-circle red line) and a variable friction coefficient p(I) (solid-cross green line), for 
different slopes 6 and erodible bed thicknesses hi. 



Figure 7: Deposit obtained in the experiments (solid-circle blue line), with the Multilayer Shallow Model using a constant 
friction coefficient p s (dotted-circle red line) and a variable friction coefficient fx(I) (solid-cross green line) and with the 
monolayer model with p s (dashed red line) and p(I) (solid green line) for a slope 0 = 22° and an erodible bed thickness 
hi = 1.82 mm. 
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Figure 8: Influence of the thickness of the erodible bed, on the runout distance rj and on the maximum final thickness 
hf (smaller graphs) with the Multilayer Shallow Model using a constant friction coefficient p, s (left hand side) and with a 
variable friction coefficient p(I) (right hand side), for a slope 9 = 22°. 


9 > 9 C , where 9 C G [12°,16°] is a critical slope, the runout distance rj and the stopping time 
tf both increase while the maximum final thickness hf decreases. Note that there is no trend 
concerning the runout when the thickness hi is increased for slopes 9 < 9 C (6 = 0°, 10°) in the 
laboratory experiments. 

Figure [T0| shows that the increase of runout distance observed in the experiments for increasing 
hi is qualitatively well reproduced with the //(/) Multilayer Shallow Model. With the //(/) Multi¬ 
layer Shallow Model, the runout increase with hi is actually larger for higher slopes, as observed 
experimentally: at 8 = 16°, the runout distance is almost unaffected by the thickness of the erodi¬ 
ble bed while it increases by 26.9% at 8 = 22° when the thickness of the erodible bed increase 
from 1.82 mm to 4.6 mm. Note that in the fi(I) MSM, the increase of the runout distance appears 
on slopes 8 > 16°, higher than 9 C in the experiments. Actually, it appears starting with the slope 
8 = 18°. When using the /%/) monolayer model, the runout distance is higher than for the Multi¬ 
layer Shallow Model whatever the slope and thickness of the erodible bed. Based on the values of 
the runout distance in these cases, it is hard to discriminate which of the monolayer or multilayer 
models is closer to the experiments. However, in the p,(I) monolayer model, the runout distance 
at 8 = 16°, 19° decreases when hi increases, contrary to the experimental data. For 8 = 22° and 
8 = 23.7°, the monolayer and multilayer //(/) models reproduce qualitatively the increase in runout 
with hi. Note that for 0 = 0°, 10° ( 8 < 8 C ), the /j(/) models predict a very slight decrease in the 
runout distance. 

As a result, the Multilayer Shallow Model with the //(/) rheology provides the results that 
are the closest to observations even though the effect of erosion is still much smaller than in the 
experiments (the runout distance increases by 4.4% for a slope 9 = 22° and from 1.82 mm to 4.6 
mm of thickness of the erodible bed, while it increases by 26.9% in the experiments). 


In Figure 11 


the final time (time at which the front stops) is plotted as a function of the 
thickness of the erodible bed for 8 = 16°, 8 = 19° and 8 = 22°. Moreover, for 8 = 22°, we also plot 
the experimental data. Experimental data show that the final time increases when the thickness of 
the erodible bed increases. In Figure [TT^l, we can see that this is true for all the values of 8 for the 
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Figure 9: Influence of the thickness hi of the erodible bed on the final runout rf for slopes 8 = 

0°, 10°, 16°, 19°, 22°, 23.7° observed in the experiments of \2$ (solid-circle blue line) and obtained with different 
simulations using the Multilayer Shallow Model with a constant friction coefficient fi s , with 20 layers (dotted-circle red 
line) and with one layer, i.e. the Savage-Hutter model / 31 / (dashed red line). There is no laboratory data for 9 = 23.7°. 
Normalisation using h 0 = 14 cm. 


multilayer method. However, in Figure ED>, for the monolayer model, we observe that it is only 
true for the highest value, 0 = 22°. At the same time, the final time decreases when the erodible 
bed increases for 0 = 16° and 6 = 19°. 


The advantage of the multilayer models is that we obtain a variable profile of the downslope 
velocity, in contrast with the constant profile of the monolayer model. It makes it possible to ob¬ 
tain a better approximation of ||Z)(u)|| (see equation (40])). As a consequence, this improves the 
approximation of the inertial number I (see equations (6F and (39)), which is a key number in the 
variable friction coefficient with n(I). 


As the main advantage of the multilayer model is the improvement of the approximation of 
||.D(tt)||, we present two approximations that can be made with the multilayer model. First, let 
us recall that a first order approximation corresponds to the definition (40). This approximation 
considers only the leading order term, i.e. ||Z)('u)|| 2=2 , , a 

in dimensionless form, we have 


Z^H 2 = 2 


= 11 Q 




Note that 


||T , ('“)II = \H>(dzU H )' +4:(d x u H ) + 2d x wd z u H + £ 2 (d x wY 


(54) 
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Figure 10: Influence of the thickness hi of the erodible bed on the final runout rj for slopes 0 = 16°, 19°, 22°, 23.7° 
observed in the experiments of \2f | / (solid-circle blue line) and obtained with different simulations using the Multilayer 
Shallow Model with the /i(I) rheology , with 20 layers (solid-cross green line) and with one layer (solid green line) and 
with the Multilayer Shallow Model with the correction of \\D(u a )\\ (dashed black line). Normalisation using ho = 14 cm. 


We can improve the approximation of ||D(m)|| at the interfaces z = z a+ i by considering the 
approximation taking into account second order terms in the previous equation. For the numerical 
tests, we consider the following approximation ||_D(w)|| at the interfaces, 

\\D{u)\\ z=Za+ ^ « ^\\Q Hja+ i\\ 2 + (d x (u H , a +i + u H)Q )) 2 . (55) 

Note that this definition corresponds to an approximation of 


\\D(u)\\ 


\! (dzUH 


2 + 4 (d x UH) 2 ■ 


at z = z Q+ i. Nevertheless, in (|54j), the term 2 d z UHd x w is not taken into account although it is 

of the same order as 4 (d x U}{) 2 . This is because when an approximation of this term is added, we 
obtain results that are very similar to those obtained when considering (55). Furthermore adding 


this term implies an additional computational cost since pre-calculated vertical velocities are re¬ 
quired. Note that (55) is a second order correction while we have developed a first order model 


that neglects other second order terms. This correction however highlights the importance of sec¬ 
ond order terms in granular collapses over erodible beds. 
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The model corresponding to the multilayer approximation with the /./(/) rheology will hereafter 
be denoted ^(J)-MSM when ||D(-u)|| is approximated by fl40|) . W hen ||.D(ti)|| is approximated by the 
correction ( |55[ ), we denote the model n(I)- C-MSM. Figure 1 10 shows that the correction of ||D(’u)|| 
corresponding to C-MSM improves the simulation of both the runout extent and the influence 
of the erodible bed. They both increase, leading to a better agreement with laboratory experiments. 




Figure 11: Influence of the thickness hi of the erodible bed on the final time tf, for slopes 6 = 16° (dashed-diamond 
red line), 0 = 19° (dashed-cross magenta line), 22° (solid-square green line) and the values observed in the experiments 
of 1 2 for slope 9 = 22° (solid-circle blue line), obtained using (a) the Multilayer Shallow Model and (b) the monolayer 
model, with the p(I) rheology with 20 layers. Normalisation using t c = ^/hg/(g cosO) and /iq = 14 cm. 


5.2.3 Flow dynamics and velocity profiles 

Figures [12] and [I3| show the time change of the granular column thickness for a slope 6 = 22° and 
an erodible bed of thickness hi = 1.82 mm for p, s and //(I), respectively, for both the monolayer 
and multilayer models. As observed for the deposit, the difference between the thickness profiles 
simulated with the multilayer and the monolayer model is stronger for n(I) than for fj, s . The p.{I)- 
MSM makes it possible to increase the maximum thickness of the flow and decrease the thickness 
of the front. This is an important result as the shape of the front may be an indicator of the flow 
rheology When a constant coefficient fj. s is used, very similar profiles are obtained with the 

Multilayer Shallow Model and monolayer model (Savage-Hutter model). As a result, the multilayer 
approach does not significantly improve the results when a constant friction coefficient is used. 
Note that during the initial instants, the simulated mass spread faster than in the experiments. 
This is partly due to the role of initial gate removal that is not taken into account here. However, 
this effect could not explain the strong difference between the simulation and experiments (see j T8] 
for more details). The hydrostatic assumption may also be responsible for this overestimation of 
the spreading velocity (see e.g. [25]). 


Figure 14 shows that the second order correction in /i(/)-C-MSM leads to simulated deposits 
that are generally closer to the experimental observations than those calculated with /r(I)-MSM. In 
particular the deposits at 9 = 19° and 6 = 22° with hi = 4.6 mm are very well reproduced (Figure 
14 3,c,d,f). However, in some cases, p.(I)-MSM gives better results than n(I)- C-MSM, for example 
for 6 = 22° with hi = 1.82 mm. This is true for the overall dynamics as illustrated in Figure [13] that 
shows the time change of the granular column thickness. We can see that with /x(/)-C-MSM, the 
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Figure 12: Thickness of the granular mass at different times in the experiments (solid-circle blue line) and with the 
Multilayer Shallow Model using a constant friction coefficient fi s and either 20 layers (dotted-circle red line) or one layer 
(Savage-Hutter model, dashed red line), for the slope 0 = 22° and erodible bed thickness hi = 1.82 mm. 


avalanche is faster and the runout is overestimated and very similar to the runout obtained with 
the l-i{I) monolayer model. As other second order terms than those included in the //(/)-C-MSM 
model are neglected, it is not easy to draw a firm conclusion on the improvement of results when 
using second order terms. 

The Multilayer approach make it possible to obtain a normal profile of the downslope velocity. 
Figures [l5| and [l6| show the normal profiles of the downslope velocity obtained at different times 
until the mass stops, for two different configurations of slopes and erodible beds. In order to obtain 
a more accurate profile, 40 layers are used in the Multilayer Shallow Model. 

The different kind of profiles observed in Figures [15] and [TTT| are in good qualitative agreement 
with typical velocity profiles of granular flows m (see also j22j and [23]). The model predicts some 
sliding at the base of the flow as shown at x = 0.095 m in Figure [15] and at x = 0.045 m in Figure 
[lO] (green profiles), in agreement with [18]. This suggests that a friction condition at the base could 
be more appropriate than the no-slip boundary condition suggested in some studies (see [10] and 
[21]). Note that the lower layers stop before the upper layers as observed experimentally. 

Let us compare the averaged velocity obtained with the monolayer model to the average of the 
velocities over all the layers in the Multilayer Shallow Model. In Figure 15, for the green profile 


(respectively red and magenta profiles), the velocity in the monolayer model is 1.01 m/s (respec¬ 
tively 0.02 and 0.14 m/s) and 0.95 m/s (respectively 0.03 and 0.05 m/s) for the averaged velocity in 
the Multilayer Shallow Model. Note that we obtain similar values for the first and second profiles. 
For the third profile, the averaged velocities strongly differ. Actually, at this position and time, 
the velocity profile corresponds to the stopping phase for the Multilayer model but not for the 
monolayer model. As a result, the velocity obtained in the Multilayer model is smaller than that 
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Figure 13: Thickness of the granular mass at different times in the experiments (solid-circle blue line), with the p(I) 
Multilayer Shallow Model using either 20 layers (solid-cross green line) or one layer (solid green line) and with the 
correction of ||_D(u a )|| (dashed black line) for the slope 0 = 22° and erodible bed thickness hi = 1.82 mm. 


obtained in the monolayer model. Figure 17 shows the normal profile of normal velocity for the 
same configuration as Figure [15] Note that the normal velocities are always negative and that their 
absolute values are greater in the upper layers. 


6 Conclusion 

In this work, we have proposed a Multilayer Shallow Model for dry granular flows that considers 
a n(I) rheology. The Multilayer approach has been applied as in [13], thus leading to a solution 
of the resulting model that is a particular weak solution of the full Navier-Stokes equations. A 
regularisation method has been used to avoid the singularity occurring when 1111(17)11 vanishes. 
A dissipative energy inequality has been proved for this model, which is an essential feature to 
guarantee that the calculated solution is physically meaningful. 

The numerical solutions of this model have been compared to the 2D analytical solutions of 
2D infinite granular layer flowing over an inclined plane proposed by |21j . The Multilayer Shallow 
Model gives an accurate approximation of this 2D analytical solution. 

By comparing the numerical results obtained with this new model to laboratory experiments, 
we have shown that the model qualitatively and sometimes quantitatively reproduces the granular 
column collapses over inclined erodible beds performed in [23]. The increase of the runout distance 
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Figure 14: Simulated deposits at different slop es 6 and erodible bed thicknesses hi with the Multilayer Shallow Model 
using the higher order correction of ||Z?(tt a )|| (55) (dashed black line) and without this correction (solid-cross green line). 
The deposits observed in the experiments is represented by solid-circle blue lines. 


with increasing thickness of the erodible bed is only reproduced when using the Multilayer Shal¬ 
low Model with the fi(I) rheology, although this increase is significantly underestimated. To our 
knowledge, this is the first time that a model has been able to reproduce this effect. The increase 
in runout distance appears for slopes 6 > 18° whereas it is observed for slopes 6 > 16° in the 
laboratory experiments. On the other hand, when using the monolayer /u(I) rheology, the increase 
of runout distance with the thickness of the erodible bed only occurs for slopes 6 > 21° . Moreover, 
in the monolayer model for 0 = 19°, the runout distance decreases as the thickness of the erodible 
bed increases, contrary to observations. As a result, when using the n(I) rheology, the multilayer 
model significantly improves the simulated deposits at different slopes over different thicknesses of 
the erodible bed compared to the monolayer model. In particular it changes the shape of the front. 
This is an important result as the shape of the front may be an indicator of the flow rheology [23, 

m- 

When considering a constant friction coefficient, the multilayer approach only slightly changes 
the results compared to the monolayer model. Even with the Multilayer model, the use of a con¬ 
stant friction coefficient does not make it possible to reproduce the increase in runout distance 
with increasing thickness of the erodible bed. The opposite effect is observed. This confirms the 
analytical results of [12] obtained for the monolayer Savage-Hutter equations. 

An important result is that this multilayer approach allows us to obtain the normal profiles of 
the downslope and normal velocities. These profiles qualitatively agree with the typical granular 
flow profiles during the developed flow and during the stopping phase H3 


39 





















Figure 15: Normal profiles of the downslope velocity obtained with the Multilayer Shallow Model (40 layers) for 0 = 22° 
and hi = 1.82 mm during granular collapse at different positions (x = 0.095, 0.495, 0.995 m). For these positions, we 
represent the velocity profiles for different times, taken every 0.15 s (blue lines). The first selected profile (green) shows a 
profile at the beginning of the flow and the second (red) and third (magenta) profiles were measured during the stopping 
stage. The final deposit is represented by the solid brown line. 



Figure 16: Normal profiles of the downslope velocity obtained with the Multilayer Shallow Model (40 layers) for 9 = 0° 
and = 1.5 mm during granular collapse at different positions (x = 0.045,0.245 m) and times taken every 0.05 s. 
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Figure 17: Normal profiles of the normal velocity obtained with the Multilayer Shallow Model (fO layers) for 9 = 22° 
and for hi = 1.82 mm during granular collapse at different positions (x = 0.095, 0.495, 0.995 m) and times taken every 


0.2 s. 


One of the differences between the multilayer and monolayer approaches is the accuracy of the 
approximation of the strain rate and consequently of the inertial number and the fi(I) friction 
coefficient. We have seen that the ^(/)-C-MSM model, which introduces a second order correction 
to improve the approximation of the strain rate, generally improves the results. The increase in 
runout distance when the thickness of the erodible bed is increased is larger and therefore closer to 
the laboratory experiments. In addition, the critical slope above which the runout increases with 
the thickness of the bed erodible is 6 > 16°, which is closer to the value observed in the experiments 
than the critical slope predicted by the model without the second order correction. This suggests 
that the extension of this shallow model up to the second order could be an important contribution. 
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